Download US Census data on natural resource dependent jobs (total # and median income) within counties and summarize across NEP boundaries:
nep_sf <- st_read(dsn = "./data-raw", layer = "NEP_Boundaries10162018", quiet = TRUE)
states <- c("AL", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", "GA",
"ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD",
"MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ",
"NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC",
"SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY",
"PR")
metrics <- load_variables(2017, "acs5", cache = TRUE) %>%
filter(grepl("Fishing|fishing", label))
totaljobs_sf <- get_acs(geography = "county", variables = c(Ag_For_Fish_Hunt_Min = "C24050_002"),
state = states, geometry = TRUE)
nr_medinc_sf <- get_acs(geography = "county", variables = c(Med_Income = "B24031_003"),
state = states, geometry = TRUE)
ustotaljobs_sf <- st_transform(totaljobs_sf, st_crs(nep_sf))
usnrmedinc_sf <- st_transform(nr_medinc_sf, st_crs(nep_sf))
nep_jobs_intersects <- st_intersects(nep_sf, ustotaljobs_sf)
nep_medinc_intersects <- st_intersects(nep_sf, nr_medinc_sf)
nep_sel_sf <- ustotaljobs_sf[unlist(nep_jobs_intersects),]
nep_sel2_sf <- nr_medinc_sf[unlist(nep_medinc_intersects),]
nep_jobs <- st_join(nep_sf, nep_sel_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_nrmedinc <- st_join(nep_sf, nep_sel2_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_jobs_sum <- nep_jobs %>%
select(NEP_NAME, estimate) %>%
group_by(NEP_NAME) %>%
summarise(jobs = sum(estimate))
nep_medinc_sum <- nep_nrmedinc %>%
select(NEP_NAME, estimate) %>%
group_by(NEP_NAME) %>%
summarise(medinc = median(estimate))
Plot US Census natural resource dependent jobs data within the NEP boundaries:
pal <- colorNumeric(palette = "viridis", domain = nep_jobs_sum$jobs, n = 10)
nep_jobs_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Jobs = ",jobs),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal(jobs)) %>%
addLegend("bottomright",
pal = pal,
values = ~ jobs,
title = "Natural Resource Dependent Jobs</br>(Ag.,Forest.,Fishing,Hunting,Mining)",
opacity = 1)
Underlying natural resource dependent jobs estimated within NEP Project Areas:
nep_jobs2 <- nep_jobs_sum %>%
select(NEP_NAME, jobs) %>%
st_set_geometry(NULL) %>%
rename("NEP" = NEP_NAME, "Total Jobs(2015)" = jobs)
nep_jobs2
Plot US Census natural resource dependent jobs, median income data within the NEP boundaries:
pal2 <- colorNumeric(palette = "viridis", domain = nep_medinc_sum$medinc, n = 10)
nep_medinc_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Median Income = ",medinc),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal2(medinc)) %>%
addLegend("bottomright",
pal = pal2,
values = ~ medinc,
title = "Natural Resource Dependent Jobs</br>(Median Annual Income)",
opacity = 1)
Underlying natural resource dependent jobs, median incomes estimated within NEP Project Areas:
nep_medinc2 <- nep_medinc_sum %>%
select(NEP_NAME, medinc) %>%
st_set_geometry(NULL) %>%
rename("NEP" = NEP_NAME, "Median Income ($, 2015)" = medinc)
nep_medinc2
Import US BEA data on natural resource dependent jobs (total county income) and summarize across NEP boundaries:
userSpecList <- list('UserID' = '88964E07-8671-4DFB-BF65-4E680DEDA641',
'Method' = 'GetData',
'datasetname' = 'Regional',
'TableName' = 'CAINC5N',
'LineCode' = '100',
'GeoFIPS' = 'COUNTY',
'Year' = '2017')
nr_income <- beaGet(userSpecList, asTable = TRUE) %>%
filter(str_detect(GeoFips, "....0")==FALSE) %>%
mutate(GEOID = as.character(GeoFips))
nr_income$GEOID <- gsub(" ", "", nr_income$GEOID, fixed = TRUE)
nr_inc_sf <- left_join(nep_sel_sf, nr_income, by = c('GEOID' = 'GEOID'))
nep_totinc <- st_join(nep_sf, nr_inc_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_inc_sum <- nep_totinc %>%
select(NEP_NAME, DataValue_2017) %>%
group_by(NEP_NAME) %>%
summarise(totinc = sum(DataValue_2017, na.rm = TRUE))
nep_inc_sum$totinc <- ifelse(nep_inc_sum$totinc==0,NA,nep_inc_sum$totinc*1000)
Plot US BEA natural resource dependent jobs total income estimates within the NEP boundaries (where available):
pal3 <- colorNumeric(palette = "viridis", domain = nep_inc_sum$totinc, n = 10)
nep_inc_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Total Income (2017) = ",totinc),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal3(totinc)) %>%
addLegend("bottomright",
pal = pal3,
values = ~ totinc,
title = "Natural Resource Dependent Jobs</br>(Total Annual Income, if reported)",
opacity = 1)
Underlying natural resource dependent jobs, total personal income estimated within NEP Project Areas:
nep_inc2 <- nep_inc_sum %>%
select(NEP_NAME, totinc) %>%
st_set_geometry(NULL) %>%
rename("NEP" = NEP_NAME, "Total Personal Income ($, 2017)" = totinc)
nep_inc2